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Abstract 

A local and distributive algorithm is proposed to find an optimal trial wave function minimizing 
the Hamiltonian expectation in a quantum system. To this end, the quantum state of the system 
is connected to the Gibbs state of a classical system with the set of couplings playing the role of 
variational parameters. The average energy is written within the replica-symmetric approximation, 
and the optimal parameters are obtained by a heuristic message-passing algorithm based on the 
Bethe approximation. The performance of this approximate algorithm depends on the structure 
and quality of the trial wave functions, starting from a classical system of isolated elements, i.e., 
mean-field approximation, and improving on that by considering the higher-order many-body in- 
teractions. The method is applied to some disordered quantum Ising models in transverse fields, 
and the results are compared with the exact ones for small systems. 
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I. INTRODUCTION 



The cavity method, relying on the Bethe approximation [1], was originally introduced 
as an alternative to the replica method to study equilibrium properties of disordered and 
effectively mean-field classical systems 

as 

. Later on it was also considered powerful 
message-passing algorithm to solve for the solutions in single instances of some computation- 
ally difficult problems {4-6]. Recently, we applied the method to a class of more challenging 
optimization problems, where the objective function itself is a computationally complex 
function, e.g., the average of a minimum energy function in a stochastic optimization prob- 
lem {7, 8]. In this study we are using these advancements to develop a message-passing 
algorithm for variational quantum-mechanics problems. 

Cavity-like approaches to quantum systems differ in nature and scope: The quantum 



belief propagation algorithms |9-|ll| are the quantum generalization of the classical belief 
propagation algorithm 12| working with local cavity density matrices instead of the cavity 
marginals. The smaller the temperature is, of course, the larger the density matrices needed 
are as quantum correlations prevail in the system. On the other hand, the quantum cav- 
ity method 14Hl8| maps the quantum problem to a classical one using the Suzuki- Trotter 



transformation, and then exploits the cavity method to estimate the relevant average quan- 
tities. This can be regarded as a dynamical mean-field theory extended to take the spatial 
correlations into account. 

In this work we take another approach that merges the techniques we used in the stochas- 



tic optimization problems 



8J with the variational principles of quantum mechanics, i.e., 
the fact that any trial wave function (density matrix) provides an upper bound for the 
ground-state energy (free energy). Alternatively, a lower bound for the free e nerg y can be 



13j. In both 



obtained by approximating the entropy with an overestimated entropy function 
cases, the problem of finding the optimal state can be recast as an optimization problem 
with an objective function that evaluates the (free) energy of a physical state. In the rest of 
this paper we shall focus on the simpler problem of finding an optimal wave function, i.e., 
at zero temperature. 

The strategy of finding the optimal wave function in a quantum system is analogous 
to that of finding the optimal configuration in a classical optimization problem. Both the 
problems are, in general, computationally hard, making efficient and accurate heuristic 



algorithms central to the study of these problems. The quantum problem is intractable 



already in one dimension 



19] or for fermionic systems due to the sign problem 20J. In 



addition, it is important for the efficiency of the variational method, to have a succinct 
representation of the trial wave functions that accurately describes the ground state of a 
quantum system 



21 



22]. 



In the following we shall map the quantum wave function to the Gibbs state of a classical 
system, considering the set of couplings as the variational parameters. The first step of 
the algorithm is to choose an appropriate classical system that best captures the quantum 
nature of the original system. The minimal set of couplings is that of a classical system 
of isolated elements, and in the maximal set one would have the whole set of many-body 
interactions. The former is equivalent to the mean-field approximation and the latter to an 
exact treatment of the problem. However, as we will see, already the two-body interactions 
give a reasonable estimate of the physical quantities in a disordered quantum Ising model. 
One may compare this with the Jastrow trial wave function, which is the product of pair 



functions 



23]- 



The second step of the algorithm is to write the objective function, which is the quantum 
average of the Hamiltonian, in terms of classical and local average quantities estimated within 
the replica-symmetric approximation ^|. The quality of this approximation depends on the 
structure of the classical interaction graph; it is expected to work well on random and sparse 
graphs, that is in effectively mean-field systems. In general, however, this approximation 
spoils the upper bound property of the Hamiltonian expectation we started from. The third 
and last step of the algorithm is to find the optimal couplings and we do this by a heuristic 
message-passing algorithm based on the Bethe approximation. 



This paper is organized as follows. We start in section [UJ with the variational quantum 
problem in its general form and write it as a classical optimization problem amenable to the 
cavity method. In section III II we write the cavity equations for a quantum spin model and 
compare the different levels of approximations with the exact results. Section [IV] gives the 
concluding remarks. 
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II. GENERAL ARGUMENTS 



Given a Hamiltonian H and a trial wave function \i/)(P)), we have (if}{P)\H\if)(P)) > E g 
where E g is the ground-state energy of H and P denotes a set of parameters characterizing 
the trial wave function. To find the optimal parameters we define the following optimization 
problem: 

z = y^e-Avtwcawcp)^ (i) 
p 

where eventually one is interested in the limit (3 opt —> oo. Note that (3 opt is just the inverse 
of a fictitious temperature and has nothing to do with the physical temperature. 

Assume H = Hq + Hi, where Hq is diagonal in the orthonormal basis \a); H and Hi 
are Hermitian operators with real eigenvalues and orthonormal eigenvectors. The trial wave 
function is represented in this representation as \i>(P)) = J2cr a (—'—)\^)- The coefficients 
a(g_;P_) are complex numbers, and \a(g_;P_)\ 2 = a(a_; P_)a*(a_; P) is a normalized probability 
distribution over a. The average energy can be written as 

= £ H^P)\ 2 [E,{a) + Ei(a)}, (2) 

where 

E (a) = (a\H \a), E 1 (a) = Re ^^^^Wxh) J . (3) 

We are going to consider fi(a;P) = \a(a;P)\ 2 as a probability measure over vari- 
ables a in a classical system and compute the average energies within the Bethe ap- 
proximation, where the classical measure is treated as if the classical interaction graph 
£ c is a tree. For example, in the Ising model the measure is approximated by /ibp = 

Hi ^Mi) ll(ij)e£- c (/%(°i> ^VM^)/^^)]); g iven tlie local marginals ^(a*) and ^(au (T,). 
The free energy in the Bethe approximation can be written in terms of the local free-energy 
shifts AFi and AF^-, corresponding to the changes in the free energy by adding spin i and 
interaction between spins i and j, respectively. To compute these free energies we need the 
cavity marginals //^(crj), i.e., the probability of having spin i in state <7j in absence of the 
interaction with spin j. The equations governing these local cavity marginals are called 



the belief propagation (BP) equations 12J. Having the cavity marginals one can obtain 



the Bethe estimation of the local marginals for any subset of the variables. Notice that in 



writing the BP equations one assumes the classical system is in a replica-symmetric (RS) 
phase. Even in a replica-symmetry-broken (RSB) phase, the approximation is still valid in 
a single pure state thanks to the exponential decay of the correlations Q]. 

The average of any local quantity like Eq and E\ can be written as a function of the BP 
cavity marginals (or messages). Therefore, the optimization problem reads 



where the indicator function I bp ensures that the messages \i satisfy the BP equations. In 
the case of multiple BP fixed points the above partition function would be concentrated 
on the one of minimum average energy for (3 opt — > oo. More accurate average energies are 
obtained, of course, by considering replica symmetry breaking and working with a probability 
distribution of the BP fixed points. 

III. QUANTUM ISING MODEL 

As an example in the following we consider the quantum Ising model with Hamiltonian 
H = H + Hi where H = —J2(ij)££ q Jij< 7 i< 7 j an d H 1 = — J2iLi hi a i ■ The interaction 
graph is defined by set S q , and iV is the size of the system. The a x,y ' z are the standard 
Pauli matrices. Here states \q) are the 2 N configurations of the of spins. In this case 

(a\H \a) = -Yj(ii)t£ q Jii i°i and {^\H\\g) = ~ Hi -o-^V^V For tne trial wave 
functions we take the Ising ansatz: 



with complex parameters P = {E>i, Kij\i = 1, . . . , N, (ij) G £ c }. By Bf and Kf- we mean 
the real part of the parameters. This results in the Gibbs measure fi(a_; P_) = \a{(r,P)\ 2 of 
a classical spin-glass model with external fields 1E>f and couplings 2K^. Notice that the 
classical interaction graph S c could be different from the quantum one S q . For simplicity, 
in the following we will assume that the two coincide as happens in zero transverse fields; 
better representations could be obtained by adding the higher order neighbors to E c . 

Given the classical measure, we write the BP equations for the cavity marginal //j_^j((7j) 




(4) 
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of spin i in the absence of spin j: 

^fa) oc e 2B ^ I I I }y 2K ^fi k ^(a k ) I = BV^, (0) 



k£di\j \ &k 

where di refers to the set of spins interacting with spin i in S c . These equations can easily 
be obtained by assuming a tree structure for the classical interaction graph 4|]. Figured] 
displays the set of variables and interactions in the classical interaction graph. Having the 
cavity marginals, the average of local energies = —J^a^a* and e$ = —hiO* read 

(eij)f, = -Jij ^2 ViVjViMh Cj), (7) 
(e,) M = -fn e~ 2B ^-^ 2K ^ cos \2Blcji + ^K^aj J ^, di {a u a di ), (8) 

with the local BP marginals //y (crj, o\,) oc e 2K ^' JlC7j fii^j(ai)jjj^i(aj) and Hifiii&i, <?di) oc 
e 2B ?°i Y[ j e9i [e 2K ^ aiai /ij-^i(aj)]. We see that the only dependence of the total average en- 
ergy on the imaginary part of the parameters comes explicitly from the (e^. Consequently, 
for hi > and hi < we can minimize the total average energy by setting B\ = 0, Kfj = 
and B\ = ir/2, K\- = 0, respectively. Without loss of generality, in the following we assume 
the hi > and the imaginary parameters are zero. 

The above average energies define the Boltzmann weight e _/3op ^- Eo ^ _/3op ^- El ^ 1 for a given 
configuration of the variational parameters and the BP messages; see figure [TJ The cavity 
marginals of the parameters (including the BP messages) can be written in a higher-level 
Bethe approximation, resembling the one-step RSB equations |4|: 

M^^,^)*^'^^ Yl e-A^" J] M k ^(K ik ,n ik )I BP (i), (9) 

Bi,{K ik ,ti ik \k£di\j} k£di\j 

where for brevity we defined //^ = /^/-n} and Jsp(i) = Yljedi^i^i-^j ~ &Pi-*j)- The 

Bethe free energy is given by F = J2i ^Fi — J2(ij)ee c ^Fij with 

e -p opt A Fi = e-^M'JlM^iK&^lBP®, (10) 

Bi^Ki^/ii^jedi} jedi 

e -p opt A Fij - e + P°^»M^(K tv ^)M^(K tv ^), (11) 

where AFi and AFy are the free-energy shifts by adding node i and link (ij), respectively. 
Notice that in the last equation we have the positive sign in the exponential to count correctly 
the energy contribution (e^)^. 



°j A 
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FIG. 1. (top) The interaction graphs of the classical system representing the trial wave function 
and (bottom) the resulting variational problem. The variables are shown with solid circles, and 
the interactions are shown with open polygons. 

The (3 opt — > oo limit of the above equations, taking the scaling Mi-+j(Kij, /%) = 



Mi-+j{Kij,iiij) = -(e i j) IJt + max {-( e i)n+ 5^ M k ^i(K ik , fi ik ) > , (12) 

Bi,{Ki k ,fJ,ik\k&di\j}:I B p(i) ,~ 

kedi\] ) 

which in short we call the MaxSum-BP equations S]- The equations can be solved by 
iteration starting from random initial messages. In each iteration we have to shift M. by 
a constant to keep max^,^. Aii->j(Kij, /%) = 0. Finally, the minimum energy is given by 
F = ~ J2(ij)e£ c ^ e ij) w ith the local energy shifts 



E g = lim^^oo 



Aej = max 

Bi,{Kij,fi i: j\jGdi}:lBp(i) 



jedi J 



Adj = max {(ey)^ + Mi^Ky, ^) + Mj ., ( /\, r //,, ) } . 
Before solving the above equations we shall consider some simpler cases. 



(13) 
(14) 



A. Zero couplings: Mean field solution 

In the zeroth order of the approximation we take = for any (ij). This is a 
mean-field approximation with a factorized measure /z(cx) = Yli^ii 17 ^^ where = 
e 2Bi<Ti /[2 cosh(25,;)]. Then using equations [7] and M we obtain the average local energies: 
i e ij)[i — — Jij tanh(2_Bj) tanh(2_Bj) and (e^ = —hi/ cosh(2£>j); therefore 

(m)\H\m)) = - E J V tanh(2B«) tanh(2^) - E i • (15) 

Here we write directly the MaxSum equations that can be used to estimate the optimal 
parameters and the minimum average energy: 

Mi^j(Bi) = -(ejp + E max{-(e ik ) fl + M k ^i(B k )} . (16) 

k£oi\j 

Then we find the optimal paramters by maximizing the local MaxSum weights: 

B* = argmax \ -(e^ + V max {-(e^ + Mj^(Bj)} \ . (17) 
Bl I jedi B] J 

B. Zero fields: Symmetric solution 

As long as the fields Bi are zero we have always a symmetric solution ^^(crj) = 1/2 
to the BP equations in the classical system. This gives the average local energies: (e^)^ = 
-Jij tanh(2ify) and (e^ = -h i /[Yl j£di cosh(2K ij )}, and 

J y taoh(2A y J - ^ n - ' 

The resulting MaxSum equations are 



wmmm) = - E ^tanh^) - e n — v ( 18 ) 



Mi^j(Kij) = -{e ij ) IM + max { -(e^ + V M k ^(K ik ) > , (19) 

and the optimal couplings are estimated by 

K*j = arg max { (e^ + (A%) + A4y .,( A',/)} . (20) 

Notice that here we have (erf) = 0, which is not the case in the ordered phase. In 
addition, the symmetric solution does not give an accurate average energy when replica 
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symmetry is broken, which may happen for large couplings in the classical system; the Bethe 
approximation works well when distant spins are nearly independent, whereas the symmetric 
solution does not respect this property in an RSB phase. To get around this problem one 
can work with the nontrivial BP fixed points, e.g., by demanding a total magnetization of 
magnitude greater than 5m C 1. At the same time one may need to limit the range of 
couplings to \Kij\ < K max in order to avoid dominance by very large couplings. 

C. General solution 

In general to solve the MaxSum-BP equations we have to work with discrete fields 
Bi G {15B\1 = —Lb, ... , Lb}, couplings Ki G {ISK\1 = —Lk,--.,Lk}, and BP cavity 
fields ui G {I5u\l = — L v , . . . , L u }. The BP cavity fields are defined by /i^fai) oc 
e Vi ^i ai . An exhaustive solution of the MaxSum-BP equations would take a time of order 
Nd(2L B )[(2L K )(2L u )] d where d is the maximum degree in £ c . Notice that given the cou- 
plings Kij and the input BP messages around spin i, one obtains for any value of 
Bf from the BP equations. The above equations can be solved more efficiently (for large 
degrees) by using a convolution function of four variables (needed to compute (e^ in the 
MaxSum-BP equations) resulting in a time complexity of order Nd(2L K )(2L u ) 6 ; see the 
Appendix. In the following instead we use a computationally easier but approximate way 
of solving the equations by restricting the domain of variables: We start by assigning a 
small number of randomly selected states Sij = {(K^, z/^), . . . , (K?, u^)} to each variable 
(Kij,Uij). Then we run the MaxSum-BP equations with these restricted search spaces to 
converge the equations and sort the states in according to their MaxSum-BP weights: 

w l3 (K l lp v\ 3 ) = (ey)„ + M^{K\ p v\.) + M^{K\ p v\ 3 ). (21) 

Next we update the search spaces by replacing the states having smaller weights with some 
other ones generated randomly but close to the best state observed during the algorithm. 
The above two steps are repeated to find better search spaces and therefore parameters. 
The algorithm performance would depend on the size of the search spaces, approaching the 
correct one for S — > oo. 

The search spaces at the beginning are chosen randomly therefore we introduce a tolerance 
5v to accept those BP messages that satisfy the BP equations within ±5v. One may start 

9 



from a large tolerance and decrease it slowly after each update of the search spaces. The 
mean-field solution mentioned before can provide a good initial point for the search spaces. 



D. Numerical results 



First, we consider a spin chain in uniform and positive transverse field but with random 
couplings. In figure [5] we compare the results obtained by the above approximations with 
the exact ones obtained by the modified Lanczos method {24)]. As expected, the mean- 
field ansatz is better than the symmetric solution in the ordered phase where the local 
magnetizations are nonzero. The reverse happens in the disordered phase and even in the 
ordered phase close to the transition point where the local magnetizations are still small. 
The two limiting behaviors are therefore displayed in the general solution. 

As another example we study the same model on a single instance of random regular 
graphs where each spin interacts with a fixed number of other randomly selected spins. 
The results displayed in figure |3] show similar behaviors observed above, except that in the 
ordered phase the symmetric solution gives a lower ground state energy than the exact 
one. As explained before, this is due to the poor estimation of the average energy by the 
symmetric solution when there are multiple BP fixed points. Moreover, due to the loops 
and small size of the system we find larger deviations from the exact data compared to the 
chain model. 

Similar qualitative behaviors are observed in ferromagnetic (J^ = 1) and ±J spin-glass 
(Jij = ±1 with equal probability) models on a random regular graph of degree d = 3. In 
the ferromagnetic case, a reasonable trial wave function is obtained by taking Bi = B and 
Kij = K. This simplification allows us to find easily the optimal parameters, and so the 
critical field h{ erTO ~ 2.29 in the thermodynamic limit. Figure H] displays the phase diagram 
of the ferromagnetic model obtained in this way. For the ± J spin-glass model, the general 
solution on single instances of size N = 1000 gives hf J ~ 2.0. The corresponding values in 



the thermodynamic limit given in Refs. 



16| and (l5| are: h{ erro ~ 2.23 and hf J ~ 1.77. 
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Chain, N=20, J=Gaussian: Lanczos — 
MF, L B =1000, 5B=0.005 




FIG. 2. Comparing the exact results with those of the mean-field (MF), the symmetric solution 
(SS) and the general solution (GS) in a single instance of the quantum Ising model in transverse 
field h on a chain of size N = 20. For comparison we also display the GS results for N = 1000. 
Here E g and m x are the ground-state energy and magnetization in the x direction, respectively. 
The inset shows the Edwards- Anderson order parameter q z = (X^(of ) 2 )/N. The couplings Jij are 
Gaussian random numbers of mean zero and variance one. SB, 5K, and 5v are the sizes of bins 
in the discrete representation of the parameters, and S is the number of states in the restricted 
domains. The data for GS are obtained by restricting the search algorithm to total magnetizations 
of magnitude greater than 6m = 0.05. 
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RRG d=3, N=20, J=Gaussian: Lanczos 
MF, L B =1000, 5B=0.005 
SS, L K =1000, 5K=0.005 




RRG d=3, N=20, J=Gaussian: Lanczos — 
MF, L B =1000, 5B=0.005 
SS, L K =1000, 5K=0.005 ■- 
GS, S=50, 5v=0.01 o 
N=1000, GS, S=30, 5v=0.02 □ 
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2 2.5 3 3.5 4 
h 



FIG. 3. The quantum Ising model in transverse field h on a random regular graph (RRG) of degree 
d = 3 and size N = 20. For comparison we also display the GS results for N = 1000. The couplings 
Jij are Gaussian random numbers of mean zero and variance one. The data for GS are obtained 
by restricting the search algorithm to total magnetizations of magnitude greater than 5m = 0.05 
and couplings of magnitude less than K max = 1. 



IV. CONCLUSION 



We suggested a heuristic message-passing algorithm to find an approximated ground state 
of a quantum system for a given instance of (possibly disordered) couplings. This was done 
by exploiting an efficient representation of the wave function and relating the quantum state 
of the system to the Gibbs state in a classical system. The local and distributive nature 
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RRG d=3: Mean Field ansatz 
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RRG d=3: Mean Field ansatz — 
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FIG. 4. The quantum Ising model with ferromagnetic interactions ( J™ = 1) in transverse field h on 
RRG of degree d = 3. The results of homogeneous mean field Bi = B and Ising {Bi = B, Ku = K) 



trial wave functions are compared with the results of the quantum cavity method 
temperature in the thermodynamic limit. 



161 ] for a small 



of the algorithm could help us to study the ground-state properties of large-scale quantum 
problems in a parallel computation. We used the following main approximations to make the 
study simple and clear: (i) working at most with the two-body interactions in the classical 
system, (ii) assuming the same classical and quantum interaction graphs, and (iii) estimating 
the classical average energies within the replica-symmetric approximation. Our approach, 
however, allows for a systematic way of improving the algorithm by considering better trial 
wave functions, using more accurate cluster variational or generalized BP approximations 
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25l . |26| . and taking the effects of replica symmetry breaking into account. Finally, the 
method can be applied to other interesting quantum systems; in a forthcoming paper we 
will show how this works in the Hubbard model in the presence of the sign problem. 
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Appendix A: Solving the MaxSum-BP equations 

Let us write the BP equations in terms of the cavity fields: 



= 2 Bl + l -Y ln[ C ° S y^ + 2 ^ ]. (Al) 
3 2/-<. cosh(z/ fc ^ - 2K ik ) 1 ; 

Using Eqs. [7] and |8] we rewrite the local average energies 

e 2K ^ coshjvj^j + Vj-tj) - e~ 2 ^' cosh(z/ i ^ j - v^) 

\ e ij)fJ.— J ij „ o K ,, „ W, . , \ . 2K,a \,/,. ,. \> v > 



13 e 2Ki 3 cosh(z/j^j + Vj-rt) + e 2K a cosh(i/j_^- - Vj^) ' 
e z ^y ii+ + e~ ZDi y ir 



o2B I7/ . , , p -2B ir ' ( A3 ) 



where 



-p-r cosh(z/ fc _^ + 2K ik ) 
1|. cosh(^) 



t7? S • ( A5 ) 



So the MaxSum-BP equations read 



e 2K zj cQshfa^j + — e 2Kij cosh(f i _ >J - — 



lJ e 2Ki i cosh(z/^j + Vj-n) + e 2Ki ^ cosh(z/,, 



+ max 



2fk 



B l ,{K lk ,v xk \kadi\j}:i BP {i) I e 2Bi y it+ + e~ 2B ^yi- 



M k ^i(K ik ,u ik ) I. (A6) 
edi\j J 
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To split the maximum over the set {K ik , fik\k £ di \j} we introduce a convolution function 
F t (x t + ,x t _,y t + ,y i L) to keep track of the quantities needed to compute the average local en- 
ergies and satisfy the BP equations. At each step t we take the maximum over the subset 
{K ikt , Vik t } for t — 1, . . . ,di — 1, where di is the degree of node % in S c , and the variables are 

1 ^ cosh(z/ A . t ,^ i + 2K ikt ,] 



+ 2^ L cosh(z/ fct/ ^-2^ fci ,) J ' 1 ; 

x* ^VM^ + ^I], (A8) 

- 2^ L cosh(z/ fct/ ^-2^ fci ,) J ' 1 ; 

_ „ cosh^^+2^) 

y+ -ll cosh^) ' (A9j 

, i-r cosh(z/fc — 2K ik ,) , 

n - ■ < Aio > 



The convolution function is updated sequentially as 

p (~t+i t+i t+i t+i\ _ 
r t +i\x + ,x_ ,y + ,y_ ) — 



max ^ F t (a;*.,a:L,S/+,yL) + ^+1^(^+1,^+1) >, (All) 



with the set of constraints I t +Y- 

4 +1 = < + J H ^w^ ], ( Ai2) 

+ + 2 L cosh(z/ fct+1 ^-2^ t+1 ) J ' 1 ; 

t+l t X i r COsh K t+i-fi + 2ififc t+1 ), OR t * 

x ^ = x In — — — ——\ = v i .- kt ,,—2B i —x\, A13 

- 2 [ cosh(z/ fct+1 ^-2^ fci+1 ) J ^ t+1 +' 1 ; 

y t+i = t cosh(u kt+1 ^ + 2K ikt+1 ) ^ 

+ + cosh(u kt+1 ^i) 

«+i = t coshK ^-2^ +1 ) 

cosh(z/ fct+1 _+i) 

and the following boundary condition: F (0,x_, 1, 1) = 0; otherwise, it is — oo. Finally, the 
MaxSum-BP message is computed after t — di — 1 steps by 

KA (K j. \ t e2 ^ J cosh(i/^j + Uj^j) - e~ 2K ^ cosh(^ - u^) 



15 e 2i ^ cosh(z/j_ > j + Vj^i) + e 2X ^ cosh(z/j. 
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